{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forecasting air pollutant concentration\n",
    "\n",
    "[Feature Engineering for Time Series Forecasting](https://www.trainindata.com/p/feature-engineering-for-forecasting)\n",
    "\n",
    "- In this notebook, we will forecast the CO concentration in the next hour. \n",
    "\n",
    "- We will train models using the features we created in the previous notebook\n",
    "\n",
    "- We will compare the performance of a naive forecast, with random forests and linear regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.metrics import mean_squared_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CO_sensor</th>\n",
       "      <th>Month</th>\n",
       "      <th>Week</th>\n",
       "      <th>Day</th>\n",
       "      <th>Day_of_week</th>\n",
       "      <th>Hour</th>\n",
       "      <th>is_weekend</th>\n",
       "      <th>CO_sensor_lag_1</th>\n",
       "      <th>RH_lag_1</th>\n",
       "      <th>CO_sensor_lag_24</th>\n",
       "      <th>RH_lag_24</th>\n",
       "      <th>CO_sensor_window</th>\n",
       "      <th>RH_window</th>\n",
       "      <th>Month_sin</th>\n",
       "      <th>Month_cos</th>\n",
       "      <th>Hour_sin</th>\n",
       "      <th>Hour_cos</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date_Time</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2004-04-05 00:00:00</th>\n",
       "      <td>1065.0</td>\n",
       "      <td>4</td>\n",
       "      <td>15</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1188.0</td>\n",
       "      <td>60.8</td>\n",
       "      <td>1224.0</td>\n",
       "      <td>56.5</td>\n",
       "      <td>1165.666667</td>\n",
       "      <td>58.566667</td>\n",
       "      <td>0.866025</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2004-04-05 01:00:00</th>\n",
       "      <td>999.0</td>\n",
       "      <td>4</td>\n",
       "      <td>15</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1065.0</td>\n",
       "      <td>65.8</td>\n",
       "      <td>1215.0</td>\n",
       "      <td>59.2</td>\n",
       "      <td>1149.666667</td>\n",
       "      <td>61.800000</td>\n",
       "      <td>0.866025</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.269797</td>\n",
       "      <td>0.962917</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2004-04-05 02:00:00</th>\n",
       "      <td>911.0</td>\n",
       "      <td>4</td>\n",
       "      <td>15</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>999.0</td>\n",
       "      <td>79.2</td>\n",
       "      <td>1115.0</td>\n",
       "      <td>62.4</td>\n",
       "      <td>1084.000000</td>\n",
       "      <td>68.600000</td>\n",
       "      <td>0.866025</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.519584</td>\n",
       "      <td>0.854419</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2004-04-05 03:00:00</th>\n",
       "      <td>873.0</td>\n",
       "      <td>4</td>\n",
       "      <td>15</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>911.0</td>\n",
       "      <td>80.0</td>\n",
       "      <td>1124.0</td>\n",
       "      <td>65.0</td>\n",
       "      <td>991.666667</td>\n",
       "      <td>75.000000</td>\n",
       "      <td>0.866025</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.730836</td>\n",
       "      <td>0.682553</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2004-04-05 04:00:00</th>\n",
       "      <td>881.0</td>\n",
       "      <td>4</td>\n",
       "      <td>15</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>873.0</td>\n",
       "      <td>81.0</td>\n",
       "      <td>1028.0</td>\n",
       "      <td>65.3</td>\n",
       "      <td>927.666667</td>\n",
       "      <td>80.066667</td>\n",
       "      <td>0.866025</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.887885</td>\n",
       "      <td>0.460065</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     CO_sensor  Month  Week  Day  Day_of_week  Hour  \\\n",
       "Date_Time                                                             \n",
       "2004-04-05 00:00:00     1065.0      4    15    5            0     0   \n",
       "2004-04-05 01:00:00      999.0      4    15    5            0     1   \n",
       "2004-04-05 02:00:00      911.0      4    15    5            0     2   \n",
       "2004-04-05 03:00:00      873.0      4    15    5            0     3   \n",
       "2004-04-05 04:00:00      881.0      4    15    5            0     4   \n",
       "\n",
       "                     is_weekend  CO_sensor_lag_1  RH_lag_1  CO_sensor_lag_24  \\\n",
       "Date_Time                                                                      \n",
       "2004-04-05 00:00:00           0           1188.0      60.8            1224.0   \n",
       "2004-04-05 01:00:00           0           1065.0      65.8            1215.0   \n",
       "2004-04-05 02:00:00           0            999.0      79.2            1115.0   \n",
       "2004-04-05 03:00:00           0            911.0      80.0            1124.0   \n",
       "2004-04-05 04:00:00           0            873.0      81.0            1028.0   \n",
       "\n",
       "                     RH_lag_24  CO_sensor_window  RH_window  Month_sin  \\\n",
       "Date_Time                                                                \n",
       "2004-04-05 00:00:00       56.5       1165.666667  58.566667   0.866025   \n",
       "2004-04-05 01:00:00       59.2       1149.666667  61.800000   0.866025   \n",
       "2004-04-05 02:00:00       62.4       1084.000000  68.600000   0.866025   \n",
       "2004-04-05 03:00:00       65.0        991.666667  75.000000   0.866025   \n",
       "2004-04-05 04:00:00       65.3        927.666667  80.066667   0.866025   \n",
       "\n",
       "                     Month_cos  Hour_sin  Hour_cos  \n",
       "Date_Time                                           \n",
       "2004-04-05 00:00:00       -0.5  0.000000  1.000000  \n",
       "2004-04-05 01:00:00       -0.5  0.269797  0.962917  \n",
       "2004-04-05 02:00:00       -0.5  0.519584  0.854419  \n",
       "2004-04-05 03:00:00       -0.5  0.730836  0.682553  \n",
       "2004-04-05 04:00:00       -0.5  0.887885  0.460065  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load data (we created this file in the previous notebook)\n",
    "\n",
    "data = pd.read_csv(\n",
    "    \"air_qual_preprocessed.csv\",\n",
    "    parse_dates=[\"Date_Time\"],\n",
    "    index_col=[\"Date_Time\"],\n",
    ")\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Separate into train and test\n",
    "\n",
    "We will train the model on a certain portion of the data and leave another part out to evaluate the model.\n",
    "\n",
    "Because these are time series, I will split based on time. We leave the most recent data as test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Timestamp('2004-04-05 00:00:00'), Timestamp('2005-04-04 14:00:00'))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Find minimum and maximum dates.\n",
    "\n",
    "data.index.min(), data.index.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((6428, 17), (494, 17))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We will keep the last month of data to test\n",
    "# the forecasting models.\n",
    "\n",
    "X_train = data[data.index <= \"2005-03-04\"]\n",
    "X_test = data[data.index > \"2005-03-04\"]\n",
    "\n",
    "X_train.shape, X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the target variable\n",
    "y_train = X_train[\"CO_sensor\"].copy()\n",
    "y_test = X_test[\"CO_sensor\"].copy()\n",
    "\n",
    "# remove raw time series from predictors set\n",
    "X_train = X_train.drop(\"CO_sensor\", axis=1)\n",
    "X_test = X_test.drop(\"CO_sensor\", axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naive Forecast\n",
    "\n",
    "- We will predict the pollutant's concentration at time t with its value at time t-1. \n",
    "\n",
    "- We will use this model as a benchmark."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train rmse:  102.99696786276506\n",
      "test rmse:  111.12015235369915\n"
     ]
    }
   ],
   "source": [
    "# Performance of naive forecast\n",
    "\n",
    "print(\"train rmse: \", mean_squared_error(\n",
    "    y_train, X_train[\"CO_sensor_lag_1\"], squared=False,))\n",
    "\n",
    "print(\"test rmse: \",  mean_squared_error(\n",
    "    y_test, X_test[\"CO_sensor_lag_1\"], squared=False,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Forests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-1 {color: black;background-color: white;}#sk-container-id-1 pre{padding: 0;}#sk-container-id-1 div.sk-toggleable {background-color: white;}#sk-container-id-1 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-1 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-1 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-1 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-1 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-1 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-1 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-1 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-1 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-1 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-1 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-1 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-1 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-1 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-1 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-1 div.sk-item {position: relative;z-index: 1;}#sk-container-id-1 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-1 div.sk-item::before, #sk-container-id-1 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-1 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-1 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-1 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-1 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-1 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-1 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-1 div.sk-label-container {text-align: center;}#sk-container-id-1 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-1 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-1\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>RandomForestRegressor(max_depth=3, n_estimators=50, random_state=0)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-1\" type=\"checkbox\" checked><label for=\"sk-estimator-id-1\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">RandomForestRegressor</label><div class=\"sk-toggleable__content\"><pre>RandomForestRegressor(max_depth=3, n_estimators=50, random_state=0)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "RandomForestRegressor(max_depth=3, n_estimators=50, random_state=0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_model = RandomForestRegressor(\n",
    "    n_estimators=50,\n",
    "    max_depth=3,\n",
    "    random_state=0,\n",
    ")\n",
    "\n",
    "rf_model.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train rmse:  99.75721333140793\n",
      "test rmse:  108.18842299455355\n"
     ]
    }
   ],
   "source": [
    "# Performance of random forests\n",
    "\n",
    "print(\"train rmse: \", mean_squared_error(\n",
    "    y_train, rf_model.predict(X_train), squared=False,))\n",
    "\n",
    "print(\"test rmse: \",  mean_squared_error(\n",
    "    y_test, rf_model.predict(X_test), squared=False,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Small uplift with respect to the Naive forecast. The predictions made by the Random Forests return smaller rmse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression\n",
    "\n",
    "Random Forest regression models are not very good at extrapolating values beyond those observed in the data. So let's now train a Lasso regression for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>#sk-container-id-2 {color: black;background-color: white;}#sk-container-id-2 pre{padding: 0;}#sk-container-id-2 div.sk-toggleable {background-color: white;}#sk-container-id-2 label.sk-toggleable__label {cursor: pointer;display: block;width: 100%;margin-bottom: 0;padding: 0.3em;box-sizing: border-box;text-align: center;}#sk-container-id-2 label.sk-toggleable__label-arrow:before {content: \"▸\";float: left;margin-right: 0.25em;color: #696969;}#sk-container-id-2 label.sk-toggleable__label-arrow:hover:before {color: black;}#sk-container-id-2 div.sk-estimator:hover label.sk-toggleable__label-arrow:before {color: black;}#sk-container-id-2 div.sk-toggleable__content {max-height: 0;max-width: 0;overflow: hidden;text-align: left;background-color: #f0f8ff;}#sk-container-id-2 div.sk-toggleable__content pre {margin: 0.2em;color: black;border-radius: 0.25em;background-color: #f0f8ff;}#sk-container-id-2 input.sk-toggleable__control:checked~div.sk-toggleable__content {max-height: 200px;max-width: 100%;overflow: auto;}#sk-container-id-2 input.sk-toggleable__control:checked~label.sk-toggleable__label-arrow:before {content: \"▾\";}#sk-container-id-2 div.sk-estimator input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 div.sk-label input.sk-toggleable__control:checked~label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 input.sk-hidden--visually {border: 0;clip: rect(1px 1px 1px 1px);clip: rect(1px, 1px, 1px, 1px);height: 1px;margin: -1px;overflow: hidden;padding: 0;position: absolute;width: 1px;}#sk-container-id-2 div.sk-estimator {font-family: monospace;background-color: #f0f8ff;border: 1px dotted black;border-radius: 0.25em;box-sizing: border-box;margin-bottom: 0.5em;}#sk-container-id-2 div.sk-estimator:hover {background-color: #d4ebff;}#sk-container-id-2 div.sk-parallel-item::after {content: \"\";width: 100%;border-bottom: 1px solid gray;flex-grow: 1;}#sk-container-id-2 div.sk-label:hover label.sk-toggleable__label {background-color: #d4ebff;}#sk-container-id-2 div.sk-serial::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: 0;}#sk-container-id-2 div.sk-serial {display: flex;flex-direction: column;align-items: center;background-color: white;padding-right: 0.2em;padding-left: 0.2em;position: relative;}#sk-container-id-2 div.sk-item {position: relative;z-index: 1;}#sk-container-id-2 div.sk-parallel {display: flex;align-items: stretch;justify-content: center;background-color: white;position: relative;}#sk-container-id-2 div.sk-item::before, #sk-container-id-2 div.sk-parallel-item::before {content: \"\";position: absolute;border-left: 1px solid gray;box-sizing: border-box;top: 0;bottom: 0;left: 50%;z-index: -1;}#sk-container-id-2 div.sk-parallel-item {display: flex;flex-direction: column;z-index: 1;position: relative;background-color: white;}#sk-container-id-2 div.sk-parallel-item:first-child::after {align-self: flex-end;width: 50%;}#sk-container-id-2 div.sk-parallel-item:last-child::after {align-self: flex-start;width: 50%;}#sk-container-id-2 div.sk-parallel-item:only-child::after {width: 0;}#sk-container-id-2 div.sk-dashed-wrapped {border: 1px dashed gray;margin: 0 0.4em 0.5em 0.4em;box-sizing: border-box;padding-bottom: 0.4em;background-color: white;}#sk-container-id-2 div.sk-label label {font-family: monospace;font-weight: bold;display: inline-block;line-height: 1.2em;}#sk-container-id-2 div.sk-label-container {text-align: center;}#sk-container-id-2 div.sk-container {/* jupyter's `normalize.less` sets `[hidden] { display: none; }` but bootstrap.min.css set `[hidden] { display: none !important; }` so we also need the `!important` here to be able to override the default hidden behavior on the sphinx rendered scikit-learn.org. See: https://github.com/scikit-learn/scikit-learn/issues/21755 */display: inline-block !important;position: relative;}#sk-container-id-2 div.sk-text-repr-fallback {display: none;}</style><div id=\"sk-container-id-2\" class=\"sk-top-container\"><div class=\"sk-text-repr-fallback\"><pre>Lasso(alpha=100, random_state=0)</pre><b>In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. <br />On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.</b></div><div class=\"sk-container\" hidden><div class=\"sk-item\"><div class=\"sk-estimator sk-toggleable\"><input class=\"sk-toggleable__control sk-hidden--visually\" id=\"sk-estimator-id-2\" type=\"checkbox\" checked><label for=\"sk-estimator-id-2\" class=\"sk-toggleable__label sk-toggleable__label-arrow\">Lasso</label><div class=\"sk-toggleable__content\"><pre>Lasso(alpha=100, random_state=0)</pre></div></div></div></div></div>"
      ],
      "text/plain": [
       "Lasso(alpha=100, random_state=0)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linear_model = Lasso(alpha=100, random_state=0)\n",
    "\n",
    "linear_model.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "train rmse:  94.23796405647266\n",
      "test rmse:  96.30704028113563\n"
     ]
    }
   ],
   "source": [
    "# Performance of linear model\n",
    "\n",
    "print(\"train rmse: \", mean_squared_error(\n",
    "    y_train, linear_model.predict(X_train), squared=False,))\n",
    "\n",
    "print(\"test rmse: \",  mean_squared_error(\n",
    "    y_test, linear_model.predict(X_test), squared=False,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear regression model provides much better forecasts compared to the benchmark and to random forests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Importance')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3gAAAHLCAYAAACarzq8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/7ElEQVR4nO3dd5htVX3/8feHpijYwo1GAbGAEQsWrFiwYKxgTUSxYEuMWKM/MTGKxlijxpqIUYwaNbYoBmwxCthBRREURSyARrEhgkjx+/tj7+GeO8ydO3dmzuyz132/nmce5ux9Zs73MDPnns9ea31XqgpJkiRJ0vhtNXQBkiRJkqTVYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZqQ5G+T/NvQdUiStBxxHzxJ0mpJ8gPg6sAlE4f3qKofr/B7Pq6q/mdl1Y1PksOA61fVQUPXIkkaB0fwJEmr7X5VtcPEx7LD3WpIss2Qj79cY61bkjQsA54kaeqSXDnJW5L8JMlZSV6UZOv+3PWS/G+SXyT5eZL/SHKV/tw7gF2BjyT5bZL/l2TfJGfO+/4/SHL3/vPDkrw/yTuT/AZ49GKPv0CthyV5Z//5bkkqycFJzkjyqyR/leRWSb6R5NdJXj/xtY9O8rkkr09yTpJvJ7nbxPlrJjkyyS+TnJbk8fMed7LuvwL+FviL/rl/vb/fwUm+leTcJKcn+cuJ77FvkjOT/E2Sn/XP9+CJ89sneWWSH/b1fTbJ9v252yb5fP+cvp5k32X8qCVJAzPgSZLWwtuAi4HrAzcH7gE8rj8X4CXANYEbArsAhwFU1SOAH7F+VPDlS3y8A4D3A1cB/mMTj78UtwF2B/4C+Gfg74C7AzcC/jzJnefd93vATsDzgQ8muVp/7j3Amf1zfTDw4iR33UjdbwFeDPxn/9z36u/zM+C+wJWAg4FXJ7nFxPe4BnBl4FrAY4E3JLlqf+6fgFsCtweuBvw/4A9JrgUcBbyoP/5M4ANJ1m3G/yNJ0gww4EmSVtuH+lGgXyf5UJKrA/cGnlZV51XVz4BXAw8FqKrTquqTVfX7qjobeBVw541/+yX5QlV9qKr+QBeENvr4S/QPVXVBVX0COA94d1X9rKrOAo6jC41zfgb8c1VdVFX/CZwK3CfJLsA+wLP773Ui8G/AIxequ6p+t1AhVXVUVX2vOscAnwDuOHGXi4AX9o9/NPBb4AZJtgIeAzy1qs6qqkuq6vNV9XvgIODoqjq6f+xPAif0/98kSSPi/H5J0mq7/2RDlCS3BrYFfpJk7vBWwBn9+asDr6ELKTv25361whrOmPj82os9/hL9dOLz3y1we4eJ22fVhh3Mfkg3YndN4JdVde68c3tvpO4FJbkX3cjgHnTP4wrASRN3+UVVXTxx+/y+vp2Ay9ONLs53beAhSe43cWxb4NObqkeSNFsMeJKkaTsD+D2w07zgMefFQAE3qapfJrk/8PqJ8/PbPZ9HF2oA6NfSzZ9KOPk1m3r81XatJJkIebsCRwI/Bq6WZMeJkLcrcNbE185/rhvcTnI54AN0o34frqqLknyIbprrpvwcuAC4HvD1eefOAN5RVY+/zFdJkkbFKZqSpKmqqp/QTSN8ZZIrJdmqb6wyNw1zR7pphOf0a8GeNe9b/BS47sTt7wCXT3KfJNsCzwUut4LHX21/DDwlybZJHkK3rvDoqjoD+DzwkiSXT3JTujVy71zke/0U2K2fXgmwHd1zPRu4uB/Nu8dSiuqnq74VeFXf7GXrJLfrQ+M7gfsl+bP++OX7hi07b/7TlyQNyYAnSVoLj6QLJ6fQTb98P/An/bkXALcAzqFr9PHBeV/7EuC5/Zq+Z1bVOcBf061fO4tuRO9MFrfY46+2L9E1ZPk58I/Ag6vqF/25A4Hd6Ebz/gt4/ib293tf/99fJPlqP/L3FOC9dM/jYXSjg0v1TLrpnMcDvwReBmzVh88D6Lp2nk03ovcsfJ8gSaPjRueSJK2SJI+m25T9DkPXIknaMnllTpIkSZIaYcCTJEmSpEY4RVOSJEmSGuEIniRJkiQ1woAnSZIkSY0Y3UbnO+20U+22225DlyFJkiRJg/jKV77y86pat9C50QW83XbbjRNOOGHoMiRJkiRpEEl+uLFzTtGUJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEdsMXYAkaeV2O/SoNXmcH7z0PmvyOJIkaXkcwZMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRkwt4CV5a5KfJfnmRs4nyWuTnJbkG0luMa1aJEmSJGlLMM0RvLcB91zk/L2A3fuPJwD/MsVaJEmSJKl5Uwt4VXUs8MtF7nIA8PbqfBG4SpI/mVY9kiRJktS6IdfgXQs4Y+L2mf2xy0jyhCQnJDnh7LPPXpPiJEmSJGlsRtFkpaoOr6q9q2rvdevWDV2OJEmSJM2kIQPeWcAuE7d37o9JkiRJkpZhyIB3JPDIvpvmbYFzquonA9YjSZIkSaO2zbS+cZJ3A/sCOyU5E3g+sC1AVf0rcDRwb+A04Hzg4GnVIkmSJElbgqkFvKo6cBPnC3jStB5fkiRJkrY0o2iyIkmSJEnaNAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNcKAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUCAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNcKAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUCAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNcKAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUCAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNcKAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUiKkGvCT3THJqktOSHLrA+V2TfDrJ15J8I8m9p1mPJEmSJLVsagEvydbAG4B7AXsCBybZc97dngu8t6puDjwUeOO06pEkSZKk1k1zBO/WwGlVdXpVXQi8Bzhg3n0KuFL/+ZWBH0+xHkmSJElq2jZT/N7XAs6YuH0mcJt59zkM+ESSJwNXBO4+xXokSZIkqWlDN1k5EHhbVe0M3Bt4R5LL1JTkCUlOSHLC2WefveZFSpIkSdIYTDPgnQXsMnF75/7YpMcC7wWoqi8Alwd2mv+Nqurwqtq7qvZet27dlMqVJEmSpHGbZsA7Htg9yXWSbEfXROXIeff5EXA3gCQ3pAt4DtFJkiRJ0jJMLeBV1cXAIcDHgW/Rdcs8OckLk+zf3+1vgMcn+TrwbuDRVVXTqkmSJEmSWjbNJitU1dHA0fOOPW/i81OAfaZZgyRJkiRtKYZusiJJkiRJWiUGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRiw54CW5dpK7959vn2TH6ZUlSZIkSdpcSwp4SR4PvB94U39oZ+BDU6pJkiRJkrQMSx3BexKwD/AbgKr6LvDH0ypKkiRJkrT5lhrwfl9VF87dSLINUNMpSZIkSZK0HEsNeMck+Vtg+yT7Ae8DPjK9siRJkiRJm2upAe9Q4GzgJOAvgaOB506rKEmSJEnS5ttmiffbHnhrVb0ZIMnW/bHzp1WYJEmSJGnzLHUE71N0gW7O9sD/rH45kiRJkqTlWmrAu3xV/XbuRv/5FaZTkiRJkiRpOZYa8M5Lcou5G0luCfxuOiVJkiRJkpZjqQHvacD7khyX5LPAfwKHbOqLktwzyalJTkty6Ebu8+dJTklycpJ3LblySZIkSdIGltRkpaqOT/KnwA36Q6dW1UWLfU3fiOUNwH7AmcDxSY6sqlMm7rM78Bxgn6r6VRI3T5ckSZKkZVpqF02AWwG79V9ziyRU1dsXuf+tgdOq6nSAJO8BDgBOmbjP44E3VNWvAKrqZ5tRjyRJkiRpwpICXpJ3ANcDTgQu6Q8XsFjAuxZwxsTtM4HbzLvPHv33/xywNXBYVX1sKTVJkiRJkja01BG8vYE9q6qm8Pi7A/sCOwPHJrlJVf168k5JngA8AWDXXXdd5RIkSZIkqQ1LbbLyTeAam/m9zwJ2mbi9c39s0pnAkVV1UVV9H/gOXeDbQFUdXlV7V9Xe69at28wyJEmSJGnLsNQRvJ2AU5J8Gfj93MGq2n+Rrzke2D3JdeiC3UOBh827z4eAA4EjkuxEN2Xz9CXWJEmSJEmasNSAd9jmfuOqujjJIcDH6dbXvbWqTk7yQuCEqjqyP3ePJKfQre17VlX9YnMfS5IkSZK09G0SjlnON6+qo4Gj5x173sTnBTyj/5AkSZIkrcCS1uAluW2S45P8NsmFSS5J8ptpFydJkiRJWrqlNll5Pd1aue8C2wOPo9vEXJIkSZI0I5Ya8Kiq04Ctq+qSqjoCuOf0ypIkSZIkba6lNlk5P8l2wIlJXg78hM0Ih5IkSZKk6VtqSHtEf99DgPPo9rd74LSKkiRJkiRtvqUGvPtX1QVV9ZuqekFVPQO47zQLkyRJkiRtnqUGvEctcOzRq1iHJEmSJGmFFl2Dl+RA4GHAdZMcOXFqR+CX0yxMkiRJkrR5NtVk5fN0DVV2Al45cfxc4BvTKkqSJEmStPkWDXhV9cMkZwIXVNUxa1STJEmSJGkZNrkGr6ouAf6Q5MprUI8kSZIkaZmWug/eb4GTknySbpsEAKrqKVOpSpIkSZK02ZYa8D7Yf0iSJEmSZtSSAl5V/XuS7YA9+kOnVtVF0ytLkiRJkrS5lhTwkuwL/DvwAyDALkkeVVXHTq0ySZIkSdJmWeoUzVcC96iqUwGS7AG8G7jltAqTJEmSJG2eTXbR7G07F+4Aquo7wLbTKUmSJEmStBxLHcE7Icm/Ae/sbz8cOGE6JUmSJEmSlmOpAe+JwJOAuW0RjgPeOJWKJEmSJEnLstQumr9P8nrgU8Af6LpoXjjVyiRJkiRJm2WpXTTvA/wr8D26LprXSfKXVfXRaRYnSZIkSVq6zemieZeqOg0gyfWAowADniRJkiTNiKV20Tx3Ltz1TgfOnUI9kiRJkqRl2pwumkcD7wUKeAhwfJIHAlTVB6dUnyRJkiRpiZYa8C4P/BS4c3/7bGB74H50gc+AJ0mSJEkDW2oXzYOnXYgkSZIkaWWW2kXzOsCTgd0mv6aq9p9OWZIkSZKkzbXUKZofAt4CfIRuHzxJkiRJ0oxZasC7oKpeO9VKJEmSJEkrstSA95okzwc+Afx+7mBVfXUqVUmSJEmSNttSA95NgEcAd2X9FM3qb0uSJEmSZsBSA95DgOtW1YXTLEaSJEmStHxbLfF+3wSuMsU6JEmSJEkrtNQRvKsA305yPBuuwXObBEmSJEmaEUsNeM+fahWSJEmSpBVbUsCrqmOmXYgkSZIkaWUWDXhJzqXrlnmZU0BV1ZWmUpUkSZIkabMtGvCqase1KkSSJEmStDJL7aIpSZIkSZpxBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhox1YCX5J5JTk1yWpJDF7nfg5JUkr2nWY8kSZIktWxqAS/J1sAbgHsBewIHJtlzgfvtCDwV+NK0apEkSZKkLcE0R/BuDZxWVadX1YXAe4ADFrjfPwAvAy6YYi2SJEmS1LxpBrxrAWdM3D6zP3apJLcAdqmqo6ZYhyRJkiRtEQZrspJkK+BVwN8s4b5PSHJCkhPOPvvs6RcnSZIkSSM0zYB3FrDLxO2d+2NzdgRuDHwmyQ+A2wJHLtRopaoOr6q9q2rvdevWTbFkSZIkSRqvaQa844Hdk1wnyXbAQ4Ej505W1TlVtVNV7VZVuwFfBPavqhOmWJMkSZIkNWtqAa+qLgYOAT4OfAt4b1WdnOSFSfaf1uNKkiRJ0pZqm2l+86o6Gjh63rHnbeS++06zFkmSJElq3WBNViRJkiRJq8uAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUCAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNcKAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUCAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNcKAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUCAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1Ihthi5AkiRpaLsdetSaPM4PXnqfNXkcSVsuR/AkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRUw14Se6Z5NQkpyU5dIHzz0hySpJvJPlUkmtPsx5JkiRJatnUAl6SrYE3APcC9gQOTLLnvLt9Ddi7qm4KvB94+bTqkSRJkqTWTXME79bAaVV1elVdCLwHOGDyDlX16ao6v7/5RWDnKdYjSZIkSU2bZsC7FnDGxO0z+2Mb81jgo1OsR5IkSZKats3QBQAkOQjYG7jzRs4/AXgCwK677rqGlUmSJEnSeExzBO8sYJeJ2zv3xzaQ5O7A3wH7V9XvF/pGVXV4Ve1dVXuvW7duKsVKkiRJ0thNM+AdD+ye5DpJtgMeChw5eYckNwfeRBfufjbFWiRJkiSpeVMLeFV1MXAI8HHgW8B7q+rkJC9Msn9/t1cAOwDvS3JikiM38u0kSZIkSZsw1TV4VXU0cPS8Y8+b+Pzu03x8SZIkSdqSTHWjc0mSJEnS2jHgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNcKAJ0mSJEmNMOBJkiRJUiMMeJIkSZLUCAOeJEmSJDXCgCdJkiRJjTDgSZIkSVIjDHiSJEmS1AgDniRJkiQ1woAnSZIkSY0w4EmSJElSIwx4kiRJktQIA54kSZIkNWKboQuQJGlLsNuhR63J4/zgpfdZk8eRJM0mR/AkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRhjwJEmSJKkRBjxJkiRJaoQBT5IkSZIaYcCTJEmSpEYY8CRJkiSpEQY8SZIkSWqEAU+SJEmSGmHAkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRkw14CW5Z5JTk5yW5NAFzl8uyX/257+UZLdp1iNJkiRJLdtmWt84ydbAG4D9gDOB45McWVWnTNztscCvqur6SR4KvAz4i2nVJEmSVsduhx61Jo/zg5feZ00eR5JaMc0RvFsDp1XV6VV1IfAe4IB59zkA+Pf+8/cDd0uSKdYkSZIkSc2aZsC7FnDGxO0z+2ML3qeqLgbOAf5oijVJkiRJUrOmNkVzNSV5AvCE/uZvk5y6Bg+7E/DzNXictdTac/L5zLbWng+095w2+/nkZVOqZHW09vMBf0Y+n7XX2t+Rz2f2tfac1ur5XHtjJ6YZ8M4Cdpm4vXN/bKH7nJlkG+DKwC/mf6OqOhw4fEp1LijJCVW191o+5rS19px8PrOttecD7T0nn8/sa+05+XxmX2vPyecz+1p7TrPwfKY5RfN4YPck10myHfBQ4Mh59zkSeFT/+YOB/62qmmJNkiRJktSsqY3gVdXFSQ4BPg5sDby1qk5O8kLghKo6EngL8I4kpwG/pAuBkiRJkqRlmOoavKo6Gjh63rHnTXx+AfCQadawAms6JXSNtPacfD6zrbXnA+09J5/P7GvtOfl8Zl9rz8nnM/tae06DP584I1KSJEmS2jDNNXiSJEmSpDVkwJMkSZKkRhjwJEmSJKkRBryGJckCxy43RC2rJcn9kjT1e5vkagscu84QtaxEkq2TfHvoOiSAJH86dA2rIckthq5BG0pyvaFr0JYnyfXm3sMl2TfJU5JcZeCytIAkV01y00FrsMnKhpLcHtiNiQ6jVfX2wQpagSRvrarHTNzeAfhwVd1twLJWJMk7gdsBH6DbemP0gSLJ54B7VdVv+tt7Au+tqhsPW9nmS/Jh4MlV9aOha1mpJCcBG32BrKpBX7xXU5LnVdULh65jNSX5UVXtOnQdm2OBMBfgw8D96P69/uraV7W6kvxvVd116DpWIskxwM50+/0eBxxbVScNW9XKJPkscAzd8/lcVZ07cEkrkmQf4DDg2nTv5wJUVV13yLpWIsmJwN5071GPpnttuFFV3XvAspYtyTrg8Vz2PfdjNvY1syzJZ4D96Z7LV4Cf0f0tPWOIeqa6TcLYJHkHcD3gROCS/nABowx4wJlJ3lhVf53kqsBRwJuHLmolquqgJFcCDgTelqSAI4B3j/gfpBcDH0lyH+AGdL9vDx+2pGW7KnByki8D580drKr9hytp2e7b//dJ/X/f0f93rD+bxTwOGF3AS/LajZ0CrrKGpayWE4AvAr+fOPZHwKvo/i0aVTBK8o35h4A95o6P9SJJVd05yXbArYB9gaOS7FBVl5mNMSKPAO4IPAh4RZLfA8dV1dOHLWvZ3gI8ne6N9iWbuO9Y/KHfY/oBwOuq6nVJvjZ0USvwYboLCv9DGz+jK1fVb5I8Dnh7VT1/gdfANWPA29DewJ7VyLBmVT0vycuT/CtwS+ClVfWBoetaqf4P6P3A9sDTgAcAz0ry2qp63aDFLUNVHZVkW+ATwI7AA6rqOwOXtVx/P3QBq6WqfgiQZL+quvnEqUOTfBU4dJjKlifJbzZ2iu5vaYwOBv6GDQPRnAPXuJbV8BDgKcDLq+qjAEm+X1V3GbasZfsB8BvgRcDv6H7XjqMbkRytJHegC0N3pLuQ8N90z2u0qur7SS4ALuw/7gLccNiqVuScub+hhlyU5EDgUaz/G9p2wHpW6gpV9eyhi1hF2yT5E+DPgb8bvJihC5gx3wSuAfxk6EJWIskDJ25+ie5N95eBSvLAqvrgMJWtXJL96d7UXZ9upOvWVfWzJFcATgFGE/CSvI71UwADXBn4HnBIEqrqKYMVt0xVdczQNUxBkuxTVZ/rb9yeca5f/jVwq6r66fwTSc5Y+3JWxfHAN6vq8/NPJDls7ctZmar6QJKPA/+Q5DF04XW0Fxyrav9+tOFw4J+q6sgkF81dPBmxz9CNDL0EOLqqLhy2nJVL8j3g58C76Ea/nlxVfxi2qhX5dJJXAB9k4gLQyKc5Hwz8FfCPfSC/DutnlozRfye5d1UdPXQhq+SFwMfppmUen+S6wHeHKsY1eECSj9D9I7ojcDO6MDT5gjCq6WVJjljkdI11fjNAkn8H3lJVxy5w7m5V9akBylqWJI9a7HxV/fta1bJakpzL+jek29FdXTyvqq40XFUrk+SWwFvpAniAXwGPGdsbhSQvAo6sqi8vcO5lY7yS2jcouqCqzh+6ltWW5OZ0UzNvVFV/PHQ9K5HkisA/0C2BuGVV7TxwSSvSN7bYB7gT3TTNPwBfqKrRzmBI8lTgDsAuwLfp1uMdW1XfG7SwZUry6QUOVwPrP7cD9uhvnlpVFw1Zz0r07xeuSPd++yLWr5Mc7fuFWWLAA5LcebHzjY5KaIYk2R7YtapOHbqW1dJ3cT0AuG1VjWo640KSXBmgqs4ZupZpSnKjqjp56DpWU5IPVNWDhq5jc/V/QzvONWAauyR7Aberqn8dupaVSnJD4M500zRvD/yoqhZ9LzEGfTO2g4FnAjtX1dYDl6Rekn2Bf6eb+hy6MP6ohS54a+0l2ZluFtk+/aHjgKdW1ZmD1GPAW2+hq9hjvbINkGQP4F+Aq1fVjfuWrftX1YsGLm3ZktyW7g/ohnQjRFsz/hGi+wH/BGxXVddJcjPghWMbOd6YJF+bt4ZtVPq21A/isp2+RteUZCmSfLWqmmrNP6bfwSR/Rteh8VNV9YOJ44+pqrcOVtgKJNl2/khDkp2q6udD1bQSSU6nG+X6LHAs8OWxT9NM8kq6EbwdgC/QvTk9rqpOH7SwzZTkoKp6Z5IFOxdW1avWuqbVkuQrwMPmLgT37/HeXVW3HLayzZPkT6vq2wt0DQbGO402ySfppjjPTZs9CHh4Ve03RD2uwdvQfsD8MHevBY6NxZuBZwFvAqiqbyR5F92C97F6PfBQ4H10TXEeyfrpCmN1GHBrunUdVNWJ/dzt0Zm3/nMrup/RBQOVs1o+DJxDt+ZmoWYerbnM/pkNGMWVzCQvpnuT/VXgb5P880TjqEPopgqPRpK70L3ZuXzfmOgJE6H1E8BYLyRcf+Tr0xbyBbrmPpdZozsyV+z/u+OgVUzHtpOzfKrqO32DtrF5BvAE4JULnBtdt+AJ66pqconU25I8bahiDHhAkicCfw1cd15L0x2ByyzeH5ErVNWXs+F+5xcPVcxqqarTkmxdVZcAR/Rtgp8zdF0rcFFVnTPv5zTWNw+T3fEupptKcsAwpayanavqnkMXsYZGEYYadT/g5n0r9MOAdyW5bt+qfozB++XAn1XVyUkeDHwyySOq6ouM8/nMuWbfJGsmpmKthqp6f5L9k9ypP3RMVX1k0KKWoarmLmi/YOhapuCEJP8GvLO//XC6rVVGpaqe0P93rN2BN+YXSQ4C3t3fPhD4xVDFGPA67wI+StcRa3Kt0LlV9cthSloVP09yPfo3bP0/sKPuEAqc3y8yPjHJy+mezxg7Gk46OcnDgK2T7E7XJn2UFxaq6uCha5iCzye5SY18I+Mt3FjCxDZVdTFAVf26n759eJL30U1JH5vt5tZz9gHiW8AHkzybcV9IOILufcND+tsH9ccGmYq1GpK8hG4myX/0h56S5HZV9bcDlrVs/fuDue05PgbcFHh6Vb1z0S+cbU+k25d1rsP2ccAbhytnZZI8BPhYVZ2b5Ll0I/r/UFVj3dvvMXRLiF5N9/r2ebr1rINwDd48SbYGrs6Ga21+NFxFy9dP8zucbgH4r4DvAwdNrusYmyTXBn5K92bn6XSdDd9YVacNWtgK9Fs8/B1wj/7Qx4EXVdXopjbO2iLj1ZDkFLptOb5PN0VzrtPXKDdp3pQkX6yq2w5dx2pKco+q+sTQdWxKkv8GXjG/sVffAfVvq2pUF7OSnADct6r+b+LYznT7xl2vqkY5jS7JiVV1s00dG5N+9tLN5qae9u+FvjbW17m5n0e/Tcd96aYFHltVew1c2rL13Wgv6Gcvzf2MLjfWLsJJvlFVN023r+SLgFcAz6uq2wxcWhMMeBOSHEK3HuqnrJ8iN/o3cv2LwlZVde7QtayGJOsAqursoWtZTUmuMNYX6jmztsh4NfQXFS6jRrqX10YWtp8D/HBu9GhskpzEZUeEzqGbvvSiqhpsmszm6LvpUlW/W+DctarqrP7zUXQ6TXJ34Oyq+vq841cGDqmqfxymspVJ8im6EbvJqVgHV9XdhqtqZfqAt+/crKV+C5LPjPX9T5Jv9s3l/g14f1V9LMnXRx7wvgjcvap+29/eAfhEVd1+2MqWZ675VT96fFJVvWtMDbHmS7eN11Or6tf97asCr6yBtiYb1dXANfA04AZVdaOqukn/McoXN4AkV0/yFroXt3OT7JnksUPXtRzpHJbk58CpwHeSnJ3keUPXtlJJbt+PEn27v71XkrFOu1hXVUdU1cX9x9uAdUMXtRJ9kNsFuGv/+fmM+7XzjcAX6Ub330zXXOF9wKlJ7rHYF86wjwJH0a1JeTjwEbpw93/A24Yra/NU1e8WCnf9ubMmbo5ic+Oq+p/54a4/fs5kuEvygbWtbMUeA/w53e/XT4AHM+BUrFXyEuBrSd7Wv1H9CjDKAN777yTfBm4JfKq/MDy6WTHzXH4u3AH0n19hwHpW6qwkbwL+Ajg6XcfqMf/betO5cAdQVb8CBgurY/4fOQ1n0F31bcXb6Kb7XbO//R26EDtGT6eb9nerqrpaVV0VuA2wT5KnD1vair0a+DP6xbj9G6I7LfoVs+sXSQ5KsnX/cRADLjJeDUmeT9dJd66Rz7asX+Q+Rj+ma+Sxd99e++bA6XTrh14+aGXLd/eqek5VndR//B1w56p6Gd32Fq0Zy5rCpRpV1+Cq+mFV7V9V66rqj6vq/mNdyjGnqt4N3Bb4IPABuv0K/3PYqpavur1Xbw/sXd0WHecz0fAryRhnlZw3OQMjyS3p1hiO1Z/TvUf9sz4YXY2u8ztw6QjYmGw1WXM/Cj5YrxObrGzodOAzSY5ioh16jXfflJ2q6r1JngPQd2a7ZOiilukRwH41sW9SVZ3eB4hP0IWk0aqqM+Z10Rzrz2mmFhmvkgfQhaCvAlTVj5OMcu1Qb4/J6X1VdUq6fYlOn/c7OCZbJ7l1VX0ZIMmt6PbIhAY6By+gtbUVo3g+fefMjdZaVU/Z2LlZtcCU7bn10tdMcs0a6Z5kAJNN8qrqPOC8idMvAz655kWtzNOA9yX5Md1FnmvQjX6NUr8k5YMTt3/Cho0AP8W4tlJ5JfCFvikWdE2YBhsFN+Bt6Ef9x3aMs2PZfOcl+SPWd9G8LeMdody2FtgUt6rOzjj3gZl0RpLbA9U/l6cC3xq4pmXppzA2sUH7hAurqpLM/R1dcVNfMONOTvIvwHv6238BnNJPj7lo41820x4HvLVfkxLgN8Dj+p/VSwatTC2Za0m/D7AnMDfC9RDglEEqWrm5vcguT7dv6dfp/oZuSvd8bzdQXdM2uqtZVXV8kj8FbtAfOrUfnQS6UcmqGltoXcyofkZV9fa+sdTcPn4PrKpLXxeSXLWftrkmbLKygP5Nwtz85tFJt7HiXJv9VwE3Bk6mWwv1kIXWRMy6JF+tqgWv5Cx2bgyS7AS8Brg73QvaJ+gW6o5mamOLV7bnJHkmsDvdFMaX0I1SvruqXjtoYcvUN/L4a7oNtQE+R7cu7wK6vTNH+boHlzbvoKrGeiFrSVrrdDq2xgp9s4s7zDUl6i/MHTfmn0mSDwLPr347mCQ3Bg6rqgcPW9l0jP19w0Jae04+n5VxBG9C/4L2Drp5wPQNPR45hm5l8+wM/DPwp3SNOz4JHEv3pvQyo2AjsVeS3yxwPHRXHsfst1X18KGLWKHJzVZfADx/qEKm4JV04fs3dFdOn0f39zRKfROPV7L+yv2kMYe7+wA3Ai4/N9W0ql44aFHLtKlOp2MOEhvx7KEL2ExXBa4EzE0B3KE/NmY3qIm9Pqvqm0luOGRB2myjGvHaAq3pz8eAt6HDgWdU1acBkuxL12VuVC1oq+qZAOk2BN+brv59geck+XVV7TlgectSVVtv+l6j9c0kP6XbM+444LNjG4Goqn+f+zzJ0yZvN+AtfZvjT8KlI/xHA6NsiZ5kd7qRyD2ZuDhSVaNqdDEpyb/SdZO7C/BvdF0NvzxoUSvzRrq1J9+ge1MwNwvjykmeWCPY02/SEraxGNXzAV5K13Hy03Q/nzvRbbE0Zt/otxSYayD1cLrfv1b9YOgCpqC1KXmtBdY1/fkY8DZ0xblwB1BVnxn5epvt6a4yXrn/+DFw0qJfoTVXVddPsitwR+A+wBv6IH6zYStbttb+kTkryRur6q/7DllH0V34Gasj6EZYX00XiA5m/B2Vb1/dhrnfqKoXJHkl3dYJY/Vj4LFzs0eS7Am8EPh/dE0JxhaIPkrXOOpd/e2H0gXyuW0s7jdMWctTVUck+ShdJ2eAZ9fEZu4jdTDwRLo14NDNUviX4cpZuX5t+25MvNetqrf3/33gQGWJSzdpP7mq/nSRu43yIuqsMOBt6PQkf8+GmzSfPmA9y5LkcLqpSucCX6Jbj/eqtVzcqaVLsjPdov07AnvRXan/7KBF6VJV9fdJXt6PEt0SeGlVjW3frknbV9WnkqRvinNYkq/QTT0dq7lW4ecnuSbd1hx/MmA9K9Vap9O7z1t7ctLcepS+E/IYbQ2cTfc+ao8ke1TVmKduX0B30WfUHannJHkHcD3gRNZ3pS7g7UPVtAZ+MHQBS1VVlyQ5NcmuG9tiZLILaiOcojmgx9CtH5pr23pcf2xsdgUuB3wXOIuu7fGvhyxIi/oRcDzw4qr6q6GLWY4k57J+5O4KE+slA1RVXWmYypYvyeQV3i8Bf0837a+SPLCqPrjwV8683yfZCvhukkPoXiN2GLimlfrvJFcBXkG3nUXRTdUcq9Y6nTa1jUWSl9H9TE4G/tAfLka8NjfJPnTTTK/NhiNeY526vTewZzXWSbCxUcmr0r3WfZmJLSyqanSduGdxRNIumo1Kd5n3RnTr725Pt4bjl8AXqqqlBhijl2Qvuo6Gd6IL598Fjqmqtwxa2BYuyRGLnK5+Xd7o9G+uvwVcBfgHuunbL6+qLw5Z12rpQ9Dlx7aOdVJrnU7737m30l1IuHQbC7qAdJ+qeu+A5W22JKcCN62q32/yziOR5NvA04GvMLEP65i6OU/q9yJ7Sr+3WhM2Nio51i7VSe680PGqOmata1kNST4MPHljI5JrzYAHJDlysfNjvJowZ2L63+2B+wJ/VFVXGbQoXUbfuOMOdNM0DwKoqmsPWpQ04+aNsl7GiEdZm9TKNhb9+ruHjC1oLybJl6rqNpu+52xL8hG60dQdgZvRzbq4NIiP/P3ct2hwVLIVSY4Fbk73Ozf4iKRTNDu3A84A3k03FWuUixzmJHkK60fuLqJbg/d5uiuoNlmZMf3GmJej+xkdB9ypXxulGZBkD7pmA1evqhsnuSmwf1W9aODSNsvEG58FjfSNz2LNOYr10+1HpdFOp81sYwGcD5yY5FNsGB5GOZLS+3SSV9D9zUw+p68OV9Ky/NPQBUzRN4FrAE2MSs5b2rEdsC1w3hiXdPT+fugCJjmCx6VzZ/cDDgRuStcl790j3P8OgCSvopvS8/mWpie0Ksm6qjp7kfOPamzbgVFJcgzwLOBNc5sxJ/lmVd142Mo2z8amw8wZ67SYpRjb31CSz7K+0+n96DudVtUoG+FsbBuLqnrsoIUtU5JHLXR8TL9j8/VbPsxXVXXXNS9mFSR5WVU9e1PHxqDlUck5/bKiA4DbVtWhQ9fTAgPePP36jQPpFuu/oKpeP3BJ2sLNdZsbuo4tVZLjq+pWSb42EfBOHPE2FotK8oGqetDQdaymsf0NJflKVd0yyUlVdZPJY0PXthz99hU3nfjvDsBHq+qOQ9emNi30Nz/3+zdUTcu1JV2cm/x3dmxmbUTSKZq9Ptjdhy7c7Qa8FvivIWuSeqOeMtyAnye5Hv0Ld5IH08gUmY0Y7TTARYztb6i1TqdNbGOR5L1V9ecb2bidkYaHg6rqnUmesdD5qnrVWte0EkmeSNeg6LpJJjdq35FuGcTozAW4jY1KAqMMePPWUG9F1/n0goHKWbGq2nHu88kRyaHqMeABSd5O12XyaLpRu28OXJI0yWH2YT0JOBz40yRnAd+nb4TTqBZ/38b2nJ5KN6XxKXSdTu8KLDgtcCRa2cZibhPw+w5axeq6Yv/fHRe913i8C/go3RrWyal+5zawr9p+wPwppvda4NhYTK6hvphuH78DhilldfWNcD6U5Pls+Hu4ZpyiCST5A+s73kz+DxntHl5qx5inLLQkyRXp1kGdO3Qt0zS26YxL4d/Q7GhkG4vHAsdW1XeHrmW1JLl8v9l5M/r+Cldnwz3jZqKF/eaYHJUEvjdxake6XgsPH6QwbWAjI5J3rqrbDVGPI3hAVW01dA3a8iR5alW9Jsk+VfW5Re662DlNWZKrAy8GrllV90qyJ3C7hvcpHNt0xqUYxd9Qa51OF9vGIsmYt7HYFXhTkt3o9o07Fjiuqk4csqgV+maSn9J1cj4O+OzIQ/ghdBu3/5QNN6Mf3TRaGh2V7Lfxeh3dVl7Q/d49tarOHK6qFZmpEUlH8KSBzDXqaHHEpCX9nldHAH9XVXsl2Qb42lzzizHpr2i/fbErvknuUVWfWMOylm1j64bmjHD9UFPNFJIcscjpqqrHrFkxU9BvSP944JnAtapq64FLWpEku9LtxboPcG/g12NtJpXkNOA2Y92ofWNaGZUESPJJuvD6jv7QQcDDq2q/4apqhyN40nC+leS7wDXnLQafmxo8xiuNLdqpqt6b5DkAVXVxkkuGLmo5quqSJNdOsl1VXbiR+4wi3PUm1w39JfCmoQpZDUsNcGPpdFpVBy/lfiPcxuK5dCFoB+BrdAHvuEGLWqF+NGUfuoC3F3Ay8NlBi1qZM4DRjkAupLFRSYB1VTV5EehtSZ42VDErNWsjkgY8aSBVdWCSawAfB0Y19WoLc16SP2J9F83bMu43DqcDn0tyJOvXHo9utAugql4w93mS+0/eblxrnU6fCowm4AEPpJuCdRRdB8MvVNXvF/+Smfcj4HjgxVX1V0MXswpOBz6T5Cg23DNudK9zE54G3KChUclfJDkIeHd/+0C6DrtjdQTdiORD+tsH9ccGGZE04EkDqqr/A/ZKsh2wR3/41Kq6aMCytKG/AY4Erpfkc8A6uo2ax+p7/cdWtNM5D8bXKXMlWnuuo1r3WVW3SHIluiv1+wGHJ/lZVd1h4NJW4ubAHYCHJTkU+C5wzIjXGv+o/9iu/2hBa6OSj6Eb8Xo13Wva54EljfrPqJkakTTgSQPr1928nW5BboBd+ilLxw5amACoqq/0P6Mb0P18Rh3A50a5+s2mqarfDluRNK7AmuTGdFMZ70zXKe8MRj5Fs6q+nmTu4s8d6UYf7gyMMuA1+jrX1KhkVf2QtmYvzdSIpAFPGt6rgHtU1akASfage4G45aBVCYAkn6WbhnUc8Lkxhzu49M3pO4Cr9bd/Djyyqk4etLBlmLfh9PUn1rK2vo51VCNeSzC25/NSus6ZrwWOH/trAkCSE4DL0Y2iHAfcqX8DPkotvc5NaGJUMsnrWLxb8FPWsJzVNFMjknbRlAaW5Bvz34gudEzDSHIduivadwRuS3fl9LiqevqghS1Tks/TdQT9dH97X7p1N7cfsq7lSHLtxc6P8Q1qa51OlyLJ66vqkKHrWC1jaYIzKcm6qjp7kfNja4TTzOvcfGMflUzyqImbLwCeP3l+TL9ns8yAJw0syVvpOmK9sz/0cGDrsbcQb0mSP6GbrnRH4C7Aj6rqnsNWtTxJvl5Ve23q2Jgl2Qo4sKr+Y+halqMfNb7rxjqdjkVr21gsVZKvVdXNh65jNY1tO58WX+fmj0oCox+VbOFvZVZHJJ2iKQ3vicCTgLkXgeOANw5Xjib161J+Ttcd6y3Ak6vqD4t/1Uw7Pcnfs+HeQ6cPWM+y9Y0ungRci64RzieBQ+ga43wdGGXAo51Op01tY7EZWrxyPrZptM28zk04HHjGvFHJNwNjHpVs4W/lhInPLzMiORRH8KQZN8bpPi1J8lS67nK7AN+mW493bFV9b9DClinJVen+EZrr+Hcs8IKq+tVwVS1Pkg8DvwK+ANwN+GO6N6JPraoTByxtRZIs+AZhzNtAtHClfqnGNtq1FGN7Tgu8zh0HHDbG17k5jY5Kjur3alNm6XXOgCfNuFl6wdiS9eseDqbb1Hjnqtp64JJWrF/vdcWq+s3QtSxHkpOq6ib951sDPwF2raoLhq1sdYx9rc2k1t7ILabF1+wWn9PYJPkv4KtsOCp5y6p6wHBVbb4k57J+5O4KwPlzp+iaY11pkMJWwSy9zjlFU5p9XoUZUJJX0l0F3oGuK9bzGHFL9CTvAv4KuIRuY+MrJXlNVb1i2MqW5dLuhVV1SZIzWwh3jXYAbMJSmuAAz16retbQ54YuYCn6ac0bVVVjbsv/GLpRyQ/2t4/rj41KVbW0/+rMcgRPmnGzdEVoS5TkwXRdM3+6kfM3GtMb7yQnVtXNkjwcuAVwKPCVMXZtTXIJ69eoBdie7mrwqK8Et9IBcP42FsBpc6cY8TYWrTTBgfYa4SQ5m25fwncDX2Le2sGqOmaIutSuWR2RdARPmn1jW9zelKp6/ybu8g66oDQW2ybZFrg/8PqquijJKK/0tTBNdiOuOBfuAKrqM0muOGRBy3TfoQuYklaa4EB7jXCuAexHt8n0w4CjgHeP6SLcfI2PSo7erI5IGvCkgSW5H3DUIp0ZW5zu05KxBfA3AT+g6zJ5bL+X3CjX4DWsiQ6AG9uHcG4bC2B0+xT2vtd/bMWGAWl0Jhv3JLn/mBv5QDdVG/gY8LEkl6P7PftMkhdU1euHrW7Zbscio5LSQpyiKQ0syTvpXsA/ALy1qr49cEnaDC1MoU2yTVVdPHQd6rTS6XRT21hU1QEDlrdiLTXBgTZeywD6YHcfunC3G93v3lur6qwh61quft3n3KjkTWlgVFLTZ8CTZkD/RuhAui6NBRxB9wJ+7qCFaZPG9qaof/PzILo3PpfO4qiqFw5VkzZuzJ1OG97GorkNp2F8r2ULSfJ24MbA0cB7quqbA5e0qiZGJV9Bd9FnrKOSmjIDnjQjkvwR8AjgacC36JoSvLaqXjdkXVpcki9W1W2HrmOpknwMOAf4Cl0nTQCq6pWDFaUNLNTpFBhdp9NWt7FopQkOtNcIJ8kfWL8ucvIN7tgbLzU1Kqnpcw2eNLAk+9ON3F0feDtw66r6WZIrAKcABrwBJdkHOLGqzktyEF1DldfMrS8aU7jr7VxV9xy6CC1qz6r6Td/p9KP0nU7prtqPSZPbWNBOExxorBFOVW01dA2rbd6o5AtaG5XUdBjwpOE9CHh1VR07ebCqzk/y2IFq0nr/AuyVZC+6tUP/RhfE7zxoVcv3+SQ3qaqThi5EG9VKp9O9ksxNLQ2wfX971KMpNNIEB5puhNOSg+hGJZ8KPCW5tMfK2P+ONEUGPGlgVfWoRc59ai1r0YIurqpKcgDdm+23jDx43wF4dJLvA79npFOxGtdEp9OGt7GYv+H0sYxww2nYdCMc4D+Gq07Q5qikps81eNLAktyWbhrmDYHtgK2B87wqNxuSHEPXdvtg4E7Az+g6AN5k0MKWqQ8Ll7GxK/maDXY6nU1jboID7TbCkbZ0BjxpYElOAB4KvA/YG3gksEdVPWfQwgRAkmvQbZh7fFUdl2RXYN+qevvApS1bkjsAu1fVEUnWATtU1feHrksdO53Otlaa4EC7jXCkLZ3DvtIMqKrTgK2r6pKqOgKwCcaMqKr/q6pXVdVx/e0fjTzcPR94NjB3AWFb4J3DVaQFfBg4ALiYbu3N3Idmw579iN396ZrgXIeuA/IYbdAIB2ilEY60RXMNnjS885NsB5yY5OV0V1C9+DKwJJ+tqjskOZeG2m0DDwBuDnwVoKp+nGTHYUvSPHY6nW2tNMGBdhvhSFs0A540vEfQBbpDgKcDu9BNz9KAquoO/X9bCz8X9k1jCmDE7d1bZqfT2dZEExxouhGOtEVzDZ40A/p1UFTV2UPXorYleSawO7Af8BK67n/vqir3W5wRSU6h2xfTTqcjYRMcSbPEgCcNJN1mNs+nG7nbiu5N3MXA62ymoGlKsh9wD7rfuY9X1ScHLkkT7HQ622yCI2nWuc5HGs7TgX2AW1XV1arqqsBtgH2SPH3Y0tSqfkrm/1bVs4A306252XbgsjShD3K7AHftPz8f/72eJTbBkTTTHMGTBpLka8B+VfXzecfXAZ+oqpsPU5laluQrwB2BqwKfBU6gW5f38EEL06X6Tqd7Azeoqj2SXBN4X1XtM3BpApJ8s6puPHQdkrQxXhGUhrPt/HAHl67Dc0RF05KqOh94IPAvVfUQ4EYD16QNPQDYn35UqKp+DLTW7GfMPp/kJkMXIUkbYxdNaTgXLvOctBJJcjvg4cBj+2N20pstdjqdbXcAHp3EJjiSZpIBTxrO5P5DkwJcfq2L0RbjqXSbnP9XVZ2c5LrApweuSRt6b5I3AVdJ8ni6TqdvHrgmrXevoQuQpMW4Bk+SdKkkr6uqJw9dx5bOTqezLckdgN2r6oh+3fQOVfX9oeuSJDDgSZImJPlqVd1i6Dq2ZP2UzAuq6pIkNwBuAHy0qi4auDRhExxJs88mK5IkzZZjgcsluRbwMeARwNsGrUiTbIIjaaYZ8CRJmi12Op1tF1Y3/ckmOJJmkgFPkjQpQxegDTqdHtUfs9Pp7JjfBOd/sAmOpBliF01J2sIk2QGgqn67wOnXrHE5uiw7nc6wqvqnvgnOb+jWRz7PJjiSZolNViRpC5Hkr4FDgSvSjdSdC7ysqt44aGHaLHY6HZZNcCTNOqdoStIWIMlzgfsC+1bVH1XV1YC7APfqz2k87NY4LJvgSJppjuBJ0hYgyanAXlV1wbzj2wNfr6o9hqlMm8utLIY19/8/yZOB7avq5UlOrKqbDV2bJIEjeJK0paj54a4/+DvgDwPUI42VTXAkzTQDniRtGc5Kcrf5B/tjPxmgHi2fnU6HZRMcSTPNKZqStAVIciPgw8Bnga/0h/emW891QFWdPFRtuqzFOp0meXRVvW3Ni9KS2ARH0tAMeJK0BUhyfeAawB6s3zT7FOBU4CdV9b2hatN6djodP9dIShqa++BJ0pbhn4HnVNVbJw8muUl/7n4D1KQJfTfT29N1Oj29P3Zd4DVJrlZVLxq0QEnSKDiCJ0lbgCTHV9WtNnLupKq6yVrXpA3Z6bQNjuBJGppNViRpy3CVRc5tv1ZFaFF2Om2DTXAkDcqAJ0lbhhOSPH7+wSSPY33TFQ3LTqcjkWSHuUY4C3jNmhYjSfM4RVOStgBJrg78F3AhG3bR3A54QFX931C1qWOn09lnExxJY2DAk6QtSJK7ADfub55cVf87ZD1az06ns22iCc4h85vgAF+yCY6kWWHAkyRpBiT5b7pOpyfNO34T4MVVZafTAdkER9JYuAZPkqTZcPX54Q6gP7bb2pejeWyCI2kUDHiSJM2Gqyxyzk6nw7MJjqRRcKNzSZJmwwlJHl9Vb548aKfTmfEU4MNJFmyCM1hVkjSPa/AkSZoBdjqdbTbBkTQWBjxJkmaInU5nk01wJI2FAU+SJGkTkhxfVbfayLmTquoma12TJC3EJiuSJEmbdpVFztkER9LMMOBJkiRt2glJHj//oE1wJM0ap2hKkiRtgk1wJI2FAU+SJGmJbIIjadYZ8CRJkiSpEa7BkyRJkqRGGPAkSZIkqREGPEmSJElqhAFPkiRJkhphwJMkSZKkRvx/5sJuVtIesTMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's now plot the feature's importance\n",
    "# according to the linear model.\n",
    "\n",
    "# Create series with feature importance.\n",
    "tmp = pd.Series(np.abs(linear_model.coef_))\n",
    "\n",
    "# Let's add the variable names.\n",
    "tmp.index = X_train.columns\n",
    "\n",
    "# Let's make a bar plot.\n",
    "tmp.plot.bar(figsize=(15, 6))\n",
    "plt.title(\"Feature importance\")\n",
    "plt.ylabel(\"Importance\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the most important feature is the concentration of CO in the previous hour. But, the humidity in the previous hour also contributes to the CO concentration, as does the time of the day (CO_sensor_season), and the concentration of CO on the previous day and in the previous 3 hrs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "fsml",
   "language": "python",
   "name": "fsml"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
